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ABSTRACT 


A 3-D compressible Navier-Stokes solver has been 
developed and applied to 3-D viscous flow over dean and 
iced wings. This method uses a third order accurate finite 
volume scheme with flux difference splitting to model the 
inviscid fluxes, and second order accurate symmetric 
differences to model the viscous terms. The effects of 
turbulence are modeled using a k-c model. In the vicinity of 
the solid walls, the k and c values are modeled using 
Gorski's algebraic model. Sample results are presented for 
surface pressure distributions, for untapered swept clean 
and iced wings made of NACA 0012 airfoil sections. The 
leading edge of these sections are modified using a 
simulated ice shape. Comparisons with the experimental 
data obtained by Prof. Bragg at the University of Illinois are 
given. 

INTRODUCTION 

The aerodynamic characteristics of lifting surfaces 
such as wings and horizontal tails may be dramatically 
altered by the accumulation of ice over the leading edge. 
Even a relatively small amount of ice accumulation may 
lead to dramatic rise in pressure drag, loss in lift, and in 
some cases, premature stall. 

Under the support of the NASA Lewis Research 
Center, a research effort has been underway at Georgia 
Tech on the effects of icing on the aerodynamic 
performance of lifting surfaces. Both 2-D and 3-D analyses 
have been developed, and calibrated [Ref. 1-3]. These 
analyses used relatively simple algorithms to solve the 
compressible Navier-Stokes equations. For example, 
standard second order accurate central differences are 
used to model the derivatives in the governing equations. A 
fourth order low-pass filter is used at every time step for 
removing any high frequency spatial oscillations. These 
analyses also used simple algebraic turbulence models 
such as the Baldwin-Lomax model to account for the 
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effects of turbulence. Finally, for the sake of computational 
stability and efficiency, these methods used an implicit 
time-marching algorithm that permit the use of very large 
time steps. Reliable convergence to steady state in less 
than 1 CPU hour of computer time on a Cray Y/MP class of 
computer, on a 100,000 point grid. Reasonable correlation 
was observed between the calculations and the measured 
data obtained by Bragg and his coworkers [Ref. 4J. 

As these analyses were applied to a number of 
ioed wing calculations, two drawbacks of these approaches 
became apparent to the present investigators. First of all, 
the turbulence model used in these analyses is a simple 
algebraic eddy viscosity model, which uses the boundary 
layer thickness 6 (or an equivalent quantity such as the 
distance from the wall, y) as the length scale of turbulent 
•ddies. The algebraic model used the boundary layer edge 
velocity U(x) (or an equivalent quantity y lwl max ) as a 
measure of the velocity scale of turbulent eddies. Such a 
simple representation of turbulence works well only for 
simple attached or mildly separated flows. In the case of 
iced wings, the flow field contains boundary layer over the 
solid surfaces, reversed flow regions, and a free shear 
layer coming off sharp corners and turns associated with 
the ice shape. The turbulent length and velocity scales for 
free shear layers are clearly different from the wall 
bounded flows. Use of a simple eddy viscosity model to 
such complex flows lead to inaccurate results in regions of 
separation. For example, the length of the separation 
bubble was often inaccurately predicted. Potapczuk has 
discussed the limitations of the Baldwin-Lomax model, and 
has recommended phenomenological fixes to this model 
[Ref. 5]. 

The second drawback of the previous analyses is 
the irrecoverable filtering out and loss of information by the 
low-pass filters. These filters tend to smear out shock 
discontinuities, contact discontinuities, and shear layers. 
The rate at which the filtering is applied is proportional to 
lul+a , where lul is the magnitude of local flow, and a is the 
speed of sound. 

The present work improves upon the existing 
numerical analysis, presented in Ref. 3, and attempts to 
rectify the above two drawbacks. It uses a flux difference 
scheme, which dynamically and adaptively reduces the 
smearing of contact discontinuities, and vorticity waves, 
and leads to crisp resolution of free shear layers that arise 
In the flow [Ref. 6]. This method also uses a sophisticated 



k-c model to account for the effects of turbulence. In the 
vicinity of solid walls, this model uses a set of 
phenomenologically derived boundary conditions for k and 
c, as described in Ref. 7. The improved analysis allows 
transition to be detected automatically by the solver as 
regions of very low turbulent kinetic energy. 

The improved solver has been applied to dean 
and load, swept wing configurations tested by Prof. Bragg 
and his ooworkers at the University of Illinois, documented 
in the Ph. D. dissertation of Khodadoust [Ref. 8]. 

The rest of this paper is organized as follows. First, 
the Roe scheme is briefly outlined. Next, the k-e 
formulation is described. Finally, several sample 
applications of the improved solver to iced wing analysis 
are presented 

MATHEMATICAL AND NUMERICAL FORMULATION 

The 3-D compressible Navier-Stokes equations 
are given in a Cartesian coordinate system as 

d£ dF dG dH _ dR dS dT 

dt dx dy dz dx dy dz 

(D 

Here q is a vector of conserved flow properties; F, 
G and H are the invisdd flux vectors; R, S and T are 
viscous terms, which include the effects of turbulence 
through the eddy viscosity concept. 

The calculations are carried out in a curvilinear 
body-fitted coordinate system (|,ri,t). In this system, the 
governing equations may be written as 

dq dF dG dH dR dS df 

dt d% di) dt d% drj dt 

( 2 ) 


where q etc. are related to their Cartesian 
counterparts by the metrics of transformation. 


The objective of the time-marching algorithm is to 
advance the flow field at cells (or nodes (i.j.k) from a time 
level W to the next time level 'n+1'. For this purpose, In the 
present work, equation (2) was discretized as follows: 
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where the operators 6 { etc. represent standard 
central differences such as 

T - F 

, f 1 mu t-in 
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The numerical flux vector T differs from the 

A 

physical flux vector F in a manner specified by Roe [Ref. 
6]. Let q|_ and qR be two estimates to the flow field vector 
q, to the left side and right side of the cell face (i+l/2,j,k). 
Then, 

T V * ^2il + lr|A|r'(^ - q,) 

where. 
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Here T is a matrix that contains the left 


eigenvectors of the matrix dF I dq, and A contains the 
eigenvalues of this matrix. These quantities must be 
computed using special "Roe" averages of qi_ and qR. 


It may be shown that the above estimates of qi_ 
and qR lead to a spatially third order accurate formulation, 
on uniform grids and on mildly stretched grids. 


Equation (3) represents a system of non-linear 
algebraic equations for the flow properties vector q at time 
level 'rul'. Furthermore, each node (i.j+1,k) is coupled to 
its nine neighbors. Thus, this system of equations are 
highly coupled. The classical Beam-Warming approximate 
factorization scheme was used solve these coupled 
nonlinear equations, in a manner similar to that described 
in Reference 1. 


The starting point is the Boussnesq assumption for 
the turbulent stresses. 
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The turbulent viscosity ^ has to modeled. This is 

achieved by modeling two physically conceivable 
quantities, the turbulent kinetic energy k, and the 
dissipation rate c. The basic high Reynolds number K-e 
model integrates the following transport equations for k and 

«. 
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The turbulent production P, is given by: 


and the turbulent viscosity ^ is related to k and e by : 

k 2 

Pi-C^p— 

The basic k-t model constants are as follows. 

Cft -.09 
Cj - 1.43 
C 2 - 1.92 
o k - 1.0 
o t - 1.3 

As in the case of the Navier-Stokes equations, the 
k-c equations were cast on a curvilinear coordinate system 
in conservation form, and solved using an ADI procedure. 
The source terms appearing in the k and c equations were 
lagged by one time step for simplicity. The turbulence 
model equations were solved separately at the end of 
every time step, after the mean flow properties had been 
computed. This resulted in the solution of a 2 x 2 system 
that may be efficiently solved. 

The starting values for k and c were obtained 
assuming the flow to be in equilibrium. That is, the 
production and dissipation terms in the k-c equations were 
set to be equal. This results in a system of simultaneous 
algebraic equations for k and e, which may be solved to 
Initialize the turbulence field. These algebraic equations 
Involve the eddy viscosity, which was computed using the 
Baldwin-Lomax model. Note that the Baldwin-Lomax model 
was used only at the start of the calculations to initialize the 


turbulence field, and is not needed during subsequent time 
steps. 


Near solid walls, the classical high Reynolds 
number k-c model breaks down. In the present work, near 
the solid wall, the k and c were assumed to have the 
following variation: 


k-Ay 3 
is constant 

Herer the constant A was evaluted using the 
computed values of k sevral nodes away from the solid 
surface. 

RESULTS AND DISCUSSION 

Navier-Stokes calculations were carried out using 
the flow solver developed in the previous section, for swept 
and unswept, untapered wing-alone geometries tested by 
Bragg et a). The airfoil sections are made of NACA 0012 
airfoil sections. The wing tested had a chord length of 15 
inches and a semispan of 37.25 inches. To study the 
effects of icing, the leading edge was modified with a 
simulated ice shape. This ice shape corresponds to that 
measured at the NASA Lewis Icing Research Tunnel, for a 
NACA 0012 airfoil, at a freestream velocity of 130 mph, 
angle of attack of 4 degrees, icing time of 4 minutes, 
volume median droplet diameter of 20 microns, liquid water 
content of 2.1. grams per cubic meter and a temperature of 
18 degrees. Bragg's measurements were done at zero 
sweep and at a 30 degree sweep. 

The calculations used a 121 x 24 x 45 grid, with 91 
points on the airfoil at each span station, 14 spanwise 
stations over the wing, and 45 points in the direction 
normal to the wing surface. Roughly a quarter of these 
points are in the boundary layer region. Thus, the boundary 
layer is only moderately resolved in these preliminary 
calculations. Figure 1 shows the internally generated body- 
fitted grid used in the present simulations. Figure 2 shows 
a typical airfoil cross section, for the iced wing case. 

Calculations have been carried for at least 16 of 
the possible 32 combinations: swept and unswept, clean 
and iced, angle of attack of 4 degrees and 8 degrees, Roe 
scheme and central difference scheme, Baldwin-Lomax 
model and k-c model. For the sake of brevity, only a small 
subset of the calculations earned out are discussed here in 
detail. 


Figure 3 shows the surface pressure distribution 
for the dean wing at 8 degrees angle of attack at 3 
spanwise locations. 42%, 56%, 72% and 89%. In these 
circulations, the Roe solver was used, and a Baldwin- 
Lomax model was employed. Because of the increased 
operations counts associated with the Roe scheme, the 
calculations were typically 30% more expensive per grid 
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point, per time step, compared to the original central 
difference scheme. The Roe solver was also stiff, requiring 
smaller time steps than the central difference scheme. 
Nevert hel e ss , calculations similar to those shown on figure 
3 could be obtained with the Roe/Baktwin-Lomax solver in 
about 1500 iterations, requiring about 90 minutes of CPU 
time on a Cray Y/MP system using a single processor. 

Rgure 4 shows the Roe/k-c solver results for the 
same case. The results are comparable to the Baldwin- 
Lomax region, except near the leading edge region, where 
the Cp distribution shows a small plateau, indicative of a 
small amount of local separation. The k-c model behaved 
wed. The strong source terms cause the k-c equation set to 
be stiff. These equations were therefore solved using a 
smaller (local) time step than the basic flow equations. 
Upper and lower limits were placed on the k and c values, 
in order to prevent the solver from diverging at the start of 
the calculations. The flow solver had no other user 
adjustable constants, and was run using the standard k-c 
constants described in the previous section. The final 
results for the k and c fields were well within these upper 
and lower limits. Thus, these limits do not affect the quality 
of the final, converged solution. A typical k-c solution 
required about 2 hours of CPU time, on a Cray Y/MP 
system, using a single processor. Work is now underway 
on reducing this time requirement. 

Figure 5 shows the calculations and 
measurements for the swept iced wing at 4 degree angle of 
attack. In this case, the flow separates over the leading 
edge ice shape, and the pressure distribution forms a 
plateau. The Roe/k-c solver is seen to perform satisfactorily 
at all the stations. It is seen that the plateau over the upper 
surface is well captured, while the plateau over the lower 
surface is not correctly captured. The solution downstream 
of the plateau region is in excellent agreement with the 
measurements. 

Figure 6 shows some preliminary results for the 
velocity field at selected chordwise locations, for the iced 
swept wing. Both the k-c results and the original algebraic 
model results are shown. The grid used in these two cases 
was identical. It appears that the k-c model does a better 
)ob of predicting the velocity profile, both near the surface 
and near the boundary layer edge. The grid used in these 
two simulations was rather coarse in the normal direction. 
Additional calculations on finer grids are needed, before 
one can conclusively say which turbulence model is 
performing better. Nevertheless, the improved agreement 
with the k-c model is encouraging. 

CONCLUDING REMARKS 

An existing compressible Navier-Stokes solver for 
dean and load wing analysis has been upgraded through 
the use of state of the art upwind difference schemes, and 
a two equation k-c model. Encouraging results have been 


obtained for swept, dean and iced wing configurations. 
There is some additional overhead assodated with the Roe 
scheme and the k-c model. Work is needed on 
improvements to the algorithm to reduce the CPU time 
requirements. Additional systematic validation of the 
improved flow solver is now underway, prior to its 
adaptation for iced wing performance analysis. 
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Fig : 1 Computational Grid (151X42X43) for iced wing with 3(K> sweep 
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Fig 2 : Simulated Glaze ice accretion 



Figure 3. Surface pressure Distribution over a Clean Swept Wing at 8 Degree Angle of Attack, 
computed using Roe Scheme with Baldwin-Lomax Model, Theory, □ Expt. 
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Figure 4. Surface Pressure Distribution over a clean Swept Wing using the Roe scheme and the k-e 

turbulence model Theory. □ Expt 
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Figure 6. Velocity Profiles at a number of 
Stations for the Iced Swept Wing at 4 degrees 
Angle of Attack, Theory, □ Expt 





3 







9 







